Artificial pancreas

ABSTRACT

A method for controlling an insulin injection device and a system for delivering insulin in a patient in need thereof. The method includes the steps of defining a time interval; receiving the level of blood glucose; computing an insulin dose to be injected in the next time interval, and transmitting the computed insulin dose to be injected to the insulin injection device.

FIELD OF INVENTION

The present invention relates to the field of instrumentation withrelation to pancreas insufficiency, especially diabetes, morespecifically type-I diabetes. Indeed, this invention proposes a newmethod and a new system implementing a novel control strategy forcompensating hyperglycemia in a fasting scenario, while ensuringpositivity of the control and no hypoglycemic episodes. The novelcontrol strategy is also dedicated to hybrid closed-loop where thepatient chooses his bolus at meal time.

BACKGROUND OF INVENTION

Insulin was discovered almost 100 years ago. Until today, it is the onlytreatment for type-1 diabetes. This treatment consists in multiple dailyinsulin injections. Basal-bolus schemes are widely used. Bolus Advisorsare designed to help patient to compute bolus doses.

Functional Insulin Therapy

Same food, same injection, at same time of the day was an option fortype-1 diabetes treatment but it was not very satisfactory. Functionalinsulin therapy is an educational program that helps patient to computeinsulin injections. It defines tools as the insulin sensitivity factor(ISF) and the carbo ratio (CR). These tools, empirically estimated fromclinical protocols are used to compute insulin boluses depending onBlood Glucose (BG) level, Blood Glucose target, carbohydrates (CHO) inthe meal and previous boluses.

The definitions of these tools are the following:

-   -   ISF, also known as correction factor (CF), is the Blood Glucose        drop caused by 1 unit of rapid-acting insulin;    -   CR, is the amount of CHO that compensates the glycemic drop        caused by 1 unit of rapid-acting insulin.

ISF and CR allow to compute meal and correction boluses:

-   -   The Correction Bolus U_(BG), depending on the patient's CF, BG        level and BG target is:

$\begin{matrix}{U_{BG} = \frac{{BG}_{level} - {BG}_{target}}{CF}} & (1)\end{matrix}$

-   -   The Meal Bolus U_(Carb) depends on the patient's CR and the        amount of carbohydrates CHO in the meal:

$\begin{matrix}{U_{Carb} = \frac{CHO}{CR}} & (2)\end{matrix}$

These tools are used in everyday life by diabetic patients to computethe insulin bolus given by:

U _(Bol) =U _(BG) +U _(Carb)   (3)

Nevertheless, patients have difficulties in computing the correctinsulin doses because:

-   -   CF might vary with the time of the day, physical activity,        stress or illness;    -   CR varies according to meal composition.

Hence, every meal turns into a stressful math problem for most type-1diabetics.

The Bolus Wizard

Nowadays, glucometers and insulin pump include a Bolus Wizard.Physicians inform these calculators with individualized values of CR, CFor Blood Glucose Target according to the time of the day. Thus diabeticpatients only have to enter the estimated amount of CHO to obtaininsulin dose recommendations.

However, one of the most common error is over-correcting a post-mealrise in Blood Glucose. It occurs when the amount of insulin that isstill active in the body is not properly taken into account. This amountis called Insulin On Board (IOB). Most Bolus Wizard include the IOB toavoid hypoglycemia. The bolus is computed as:

U _(Bol) =U _(BG) +U _(Carb)−IOB   (4)

IOB is a function of the Duration of Insulin Action (DIA) and the amountof previous boluses. IOB is computed in different ways according to thedifferent Bolus Wizards. Nonetheless, incorrect estimation of DIAinduces mismatch in the IOB and insulin injection. As a consequence,hypoglycemia occurs when DIA is underestimated while overestimation ofDIA leads to hyperglycemia. Determination of individualized DIA remainsa critical point.

Since the past 50 years, closed-loop control of Blood Glucose in type-1diabetes, the so-called artificial pancreas (AP), remains a challenge.In 1977, the Biostator was the first realization of an artificialpancreas. Many families of controllers were designed, among which arethe Proportional-Integrate-Derivative (PID), PID with insulin feedback,Biohormonals, sliding modes, fuzzy-logic and model-predictivecontrollers (MPC). The latter became popular because they includedconstraints on the control and safety algorithms. Nowadays, closed-loopclinical trials are conducted for inpatients and outpatients.

Still, ambulatory artificial pancreas systems are not available becausemany improvements are needed. Among them, for MPC control algorithms:

-   -   the prediction horizon has to be extended;    -   the accuracy of predictions given by the model has to be        improved;    -   the individualization of the controller requires an engineering        expert work;

Consequently, there is still a need to provide an artificial pancreasthat could guarantee positivity of the control while avoidinghypoglycemic episodes.

SUMMARY

According to a first aspect, the present invention relates to a methodfor delivering insulin in a patient in need thereof.

The method comprising:

-   -   determining a time interval;    -   measuring the level of blood glucose of the patient at the final        endpoint of each time interval;    -   using a processor for computing a global insulin injection rate        to be injected at the final endpoint of each time interval, said        computing of the injection rate taking into account patient's        own parameters such as for example blood glucose and insulin on        board; and    -   delivering said computed global insulin injection rate to be        injected at the final endpoint of each time interval to the        patient.

According to one embodiment, said time interval ranges from 1millisecond to 3 hours, from 0.1 second to 1 hour or from 1 second to 15minutes.

According to one embodiment, the computed global insulin injection ratecomprises a constant insulin injection rate such as basal rate and avariable insulin injection rate.

According to one embodiment, the global insulin injection rate to beinjected at the final endpoint of each time interval is computed fortuning the velocity of decrease of glycemia of the patient withoutreaching hypoglycemia levels.

According to one embodiment of the method,

-   -   at least one window of time comprising a plurality of time        interval is defined;    -   for a window of time, a total amount of insulin to be injected        is determined; and    -   the global insulin injection rate at the final endpoint of each        time interval is function of the parameters of blood glucose and        insulin on board, as computed by the processor.

According to one embodiment, said window of time is a period of timehigher than 12 hours, can be 24-72 hours and can also last severalmonths to several years.

According to a second aspect, the invention relates to a computerprogram product, comprising a non-transitory tangible computer readablemedium having a computer readable program code embodied therein, whichis adapted to be executed to implement a method for delivering insulin.The method comprising:

-   -   defining a time interval;    -   measuring the level of blood glucose of the patient at the final        endpoint of each time interval;    -   using a processor for computing a global insulin injection rate        to be injected at the final endpoint of each time interval, said        computing of the injection rate taking into account the        parameters of blood glucose and insulin on board; and    -   delivering said computed global insulin injection rate to be        injected at the final endpoint of each time interval to the        patient.

According to a third aspect, the invention further relates to a systemfor delivering insulin, the system comprising a computer program productaccording to the second aspect of the present invention; an insulinpump; and a means for measuring the level of blood glucose of thepatient in a patient body such as a glucose sensor or continuous glucosemeasurement; wherein the system is capable to execute the methodaccording to the first aspect of the present invention. According to oneembodiment, said means for continuous glucose measurement is connectedto said computer program product.

According to a fourth aspect, the present invention also relates to acomputer-implemented method for controlling an insulin injection deviceof a diabetic user, comprising iteratively the steps of:

-   -   determining a time interval T_(s);    -   receiving the blood glucose level;    -   computing an insulin dose to be injected u(nT_(s)) in the next        time interval; and    -   transmitting the computed insulin dose to be injected to an        insulin injection device of a diabetic user.

The computed insulin dose to be injected u(nT_(s)) is computed accordingto the formula:

u(nT _(s))=k _(d) ×ũ(nT _(s))+T _(s) ×U _(Bas).

In this formula, U_(Bas) is the diabetic user's specific basal insulininjection rate and T_(s)×U_(Bas) corresponds to the basal dose which isalways injected in order to get closer to a real pancreas.

By “iteratively”, it had to understood that the steps of receiving,computing and transmitting (and optionally the step of determining) arecontinuously repeated.

ũ(nT_(s)) is a correction insulin dose which correspond to the insulindose to be injected in order to reach to a blood glucose level ofreference. This correction insulin dose is computed asũ(nT_(s))=u_(BG)(nT_(s))−IOB(nT_(s)) wherein u_(BG)(nT_(s)) is theinsulin dose needed to reach the blood glucose level x₁(nT_(s)) to ablood glucose level target x₁ _(ref) without considering previousinsulin injections; and IOB(nT_(s)) is the insulin dose still active inthe body of the diabetic user.

In this way, the computed insulin dose to be injected u(nT_(s)) isalways positive and does not need to be set at zero if the blood glucoselevel.

The k_(d) coefficient is a tuning parameter strictly positive andinferior or equal to 1. When k_(d) is strictly inferior to 1, the methodinjects only a part of the needed dose in order to spread in the timethe insulin dose to be injected. This tuning parameter acts like asafety parameter. Indeed, if the blood glucose level becomes lower thanexpected because of an error in parameters or because of a physicalactivity of the diabetic user, the shifting in time of a part of theinsulin injection permits to avoid hypoglycemia to the diabetic user.

According to one embodiment, the insulin dose needed to reach the bloodglucose level x₁(nT_(s)) to a blood glucose level target x₁ _(ref)without considering previous insulin injections is computed as:

${u_{BG}\left( {nT}_{s} \right)} = \frac{{x_{1}\left( {nT}_{s} \right)} - x_{1_{ref}}}{\theta_{2}}$

wherein x₁ _(ref) is the target blood glucose level; and θ₂ is thediabetic user's specific insulin sensitivity factor.

According to one embodiment, IOB(T_(s)) is computed as:

IOB(nT _(s))=θ₃×(x ₂(nT _(s))+x ₃(nT _(s)))

-   -   wherein        -   x₂(nT_(s)) is the plasma insulin rate;        -   x₃(nT_(s)) is the subcutaneous insulin rate; and        -   θ₃ is the diabetic user's specific insulin response time.

According to one embodiment, the insulin dose to be injected u(nT_(s))comprises a proportional component to the blood glucose levelx₁(nT_(s)), a derivative component to the blood glucose level x₁(nT_(s))and a second derivative component to the blood glucose level x₁(nT_(s)).

According to one embodiment, x₂(nT_(s)) is computed asx₂(nT_(s))=U_(Bas)−1/θ₂×{dot over (x)}₁(nT_(s)). According to oneembodiment, x₃(nT_(s)) is computed as x₃(nT_(s))=x₂(nT_(s))−θ₃×({umlautover (x)}₁(nT_(s)))/θ₂. In said embodiments, {dot over (x)}₁(nT_(s)) isthe time derivative of the blood glucose level x₁(nT_(s)) and {umlautover (x)}₁(nT_(s)) is the second time derivative of the blood glucoselevel x₁(nT_(s)).

According to an alternative embodiment, x₂(nT_(s)) and x₃(nT_(s)) aredetermined by an observer. Such observer can be an algorithm or a devicewhich measures the amount of insulin injected since a predeterminedtime.

According to one embodiment, the parameter k_(d) is strictly positiveand strictly inferior to 1. According to one embodiment, the parameterk_(d) is strictly positive and inferior or equal to 0.99, 0.95, 0.90,0.85 or 0.80 . . .

According to one embodiment, x₁ _(ref) ranges from 70 mg/L to 140 mg/L.According to one embodiment, the time interval T_(s) ranges from 1millisecond to 3 hours, from 0.1 second to 1 hour or from 1 second to 15minutes.

According to one embodiment, the method further comprises the step ofcomputing a second insulin dose u_(Carb) to be injected when an actuatoris activated, the second insulin dose u_(Carb) corresponding to the doseof insulin to be injected compensating a meal. According to oneembodiment, the actuator is activated before, during or after a meal orwhen a meal is detected. According to another embodiment, the actuatoris activated manually by the diabetic user. The latter corresponds tothe so-called hybrid closed-loop.

According to a fifth aspect, the present invention further relates to asystem for delivering insulin. The system comprises:

-   -   a processor comprising instructions to operate the        computer-implemented method according to the fourth aspect of        the present invention;    -   an insulin injection device; and    -   a sensor for measuring the blood glucose level of a diabetic        user.

According to one embodiment, said sensor is connected to the processorin order to provide to said processor the blood glucose levelx₁(nT_(s)).

According to one embodiment, the processor comprising a processor deviceand at least one memory element associated with the processor, the atleast one memory element storing processor-executable instructions that,when executed by the processor, perform a method of controlling deliveryof insulin from insulin injection device to the body of the diabeticuser according to the fourth aspect of the present invention.

The insulin injection device is controlled by the processor and is ableto inject into the patient body the insulin rate during a time intervalor the insulin dose at the end of each time interval computed by theprocessor with the method according to the fourth aspect of the presentinvention.

According to one embodiment, the insulin injection device comprising aninsulin reservoir for insulin to be delivered from the insulin injectiondevice to a body of a user.

In a sixth aspect, the invention relates to a closed-loop insulininfusion system comprising: a continuous glucose sensor that generatessensor data indicative of sensor glucose values for a user; and aninsulin infusion device to receive the sensor data generated by thecontinuous glucose sensor, the insulin infusion device comprising: aninsulin reservoir for insulin to be delivered from the insulin infusiondevice to a body of a user; a processor architecture comprising at leastone processor device; and at least one memory element associated withthe processor architecture, the at least one memory element storingprocessor-executable instructions that, when executed by the processorarchitecture, perform a method of controlling closed-loop delivery ofinsulin from the insulin reservoir to the body of the user, the methodcomprising:

-   -   initiating a closed-loop operating mode of the insulin infusion        device; in response to initiating the closed-loop operating        mode, obtaining a most recent sensor glucose value for the user;    -   computing a current insulin on board (IOB(nT_(s))) value that        indicates an amount of active insulin in the body of the user;    -   determining an insulin dose to be injected u(nT_(s)) during a        predetermined time interval (Ts);    -   operating the insulin infusion device in a closed-loop mode to        deliver insulin from the insulin reservoir to the body of the        user in accordance with insulin dose to be injected that is        determined, wherein the insulin dose to be injected represents        an amount of insulin to be delivered during each time interval        (Ts).

According to a seventh aspect, the invention relates to a method fordelivering insulin in a patient in need thereof, the method comprisingthe steps of:

-   -   defining a time interval;    -   measuring the level of blood glucose x₁(t) of the patient at the        final endpoint of each time interval;    -   using a processor for computing a global insulin injection rate        u_(i)(t) to be injected during the next time interval or at the        final endpoint of each time interval; and    -   delivering said computed global insulin injection rate u_(i)(t)        to be injected during the next time interval or at the final        endpoint of the next time interval to the patient.

In said aspect of the invention, u_(i)(t)=k×ũ₁(t)+U_(Bas).

U_(bas) is a constant patient's specific basal insulin rate; k is atuning parameter strictly positive and inferior or equal to 1; ũ₁(t) isa variable insulin injection rate computed as:

${{{\overset{\sim}{u}}_{1}(t)} = {\frac{{x_{1}(t)} - x_{1{ref}}}{\theta_{2}} - {\theta_{3}{x_{2}(t)}} - {\theta_{3}{x_{3}(t)}}}};$

wherein

-   -   x₂(t) is computed as x₂(t)=U_(Bas)−1/θ₂×{dot over (x)}₁(t);    -   x₃(t) is computed as x₃(t)=x₂(t)−θ₃×({umlaut over (x)}₁(t))/θ₂;    -   x_(1ref) is a blood glucose level target;    -   {dot over (x)}₁(t) is the time derivative of the blood glucose        level x₁(t);    -   {umlaut over (x)}₁(t) is the second time derivative of the blood        glucose level x₁(t);    -   θ₂ is the patient's specific insulin sensitivity factor;    -   θ₃ is the diabetic user's specific insulin response time.

According to one embodiment, the steps of measuring the level of bloodglucose, using a processor for computing global insulin injection rateand delivering said computed global injection rate are continuouslyexecuted at each time interval, optionally during a predetermined windowof time.

According to one embodiment, the parameter k is strictly positive andstrictly inferior to 1. According to another embodiment, the parameter kis equal to 1.

The insulin dose delivered at each time interval (according to thefourth aspect) equals the insulin rate (according to the eighth aspect)times the time interval. Thus k_(d) defines as:

-   -   k_(d)=1−e^((−k·T) ^(s) ⁾ where k is in rad/s and k_(d) is        dimensionless

According to an eighth aspect, the present invention relates to acomputer program product comprising instructions which, when the programis executed by a computer, causes the computer to carry out the stepsof:

-   -   receiving a level of blood glucose x₁(t) of a patient;    -   computing a global insulin injection rate u_(i)(t) to be        injected;    -   transmitting the computed global insulin injection rate u_(i)(t)        to an insulin injection device;

wherein u_(i)(t)=k×ũ₁(t)+U_(Bas); and U_(Bas) is a constant patient'sspecific basal insulin rate, k is a tuning parameter strictly positiveand inferior or equal to 1 and ũ₁(t) is a variable insulin injectionrate computed as:

${{{\overset{\sim}{u}}_{1}(t)} = {\frac{{x_{1}(t)} - x_{1{ref}}}{\theta_{2}} - {\theta_{3}{x_{2}(t)}} - {\theta_{3}{x_{3}(t)}}}};$

and further wherein:

-   -   x₂(t) is computed as x₂(t)=U_(Bas)−1/θ₂×{dot over (x)}₁(t);    -   x₃(t) is computed as x₃(t)=x₂(t)−θ₃×({umlaut over (x)}₁(t))/θ₂;    -   x_(1ref) is a blood glucose level target;    -   {dot over (x)}₁(t) is the time derivative of the blood glucose        level x₁(t);    -   {umlaut over (x)}₁(t) is the second time derivative of the blood        glucose level x₁(t);    -   θ₂ is the patient's specific insulin sensitivity factor;    -   θ₃ is the diabetic user's specific insulin response time.

According to one embodiment, the parameter k is strictly positive andstrictly inferior to 1. According to another embodiment, the parameter kis equal to 1.

According to a ninth aspect, the invention also relates to a system fordelivering insulin, the system comprising: a computer program productaccording to the invention; an insulin pump; and a means for measuringthe level of blood glucose of the patient in a patient body such as aglucose sensor or continuous glucose measurement; wherein the system iscapable to execute the method according to the present invention.

According to a tenth aspect, the invention relates to a method forcontrolling an insulin injection device of an user, comprisingiteratively the steps of:

-   -   determining a time interval T_(s);    -   receiving the blood glucose level x₁(nT_(s));    -   computing an insulin dose to be injected u(nT_(s)) in the next        time interval;    -   wherein the insulin dose to be injected u(nT_(s)) comprises at        least:    -   a first term being a function of the comparison between the        received blood glucose level x₁(nT_(s)) and a predefined blood        glucose level target x₁ _(ref) in a preliminary step of the        method;    -   a second term, the second term being an estimated value of the        insulin dose still active in the body IOB(nT_(s)) of the user;    -   the second term being superior or equal to the first term at        each iteration of the method.

The advantage of the feature “the second term being superior or equal tothe first term at each iteration of the method” is to preserve thepositivity of the computed insulin dose. In this way, the computedinsulin dose to be injected is always positive and does not need to beset at zero in case of an insulin dose to be injected become negative.The positivity of the command ensure a security to the user.

In one embodiment, the first and the second terms of the insulin dose tobe injected is a function of a corrective factor inferior or equal to 1,the corrective factor being configured to adapt the duration of theinjection to a predefined duration reference.

In one embodiment, the first and the second terms of the insulin dose tobe injected is a linear function of the corrective factor.

The advantage of this embodiment is to inject only a fraction of thecalculated insulin dose which the user theoretically needs to reach theblood glucose level target. This tuning parameter acts like a safetyparameter. Indeed, if the blood glucose level becomes lower thanexpected because of an error in parameters or because of a physicalactivity of the diabetic user, the shifting in time of a part of theinsulin injection permits to avoid hypoglycemia to the diabetic user.

In one embodiment, the insulin dose to be injected comprises a thirdterm calculated on at least one specific injection rate of a predefineduser profile. In one embodiment, said third term is constant on eachiteration. In one embodiment, said third term is constant along apredefined time comprising a plurality of adjacent iterations.

The advantage of said third term is to provide a basal rate of insulinto mimic the behavior of a health human pancreas and to ensure that atleast a minimum amount of insulin is injected at each iteration.

In one embodiment, the insulin dose to be injected is a function of atleast one of the following predefined user profile parameters: aspecific insulin response time; and/or a specific insulin sensitivityfactor. The advantage of said embodiment is to use some coefficientwhich is usually handle by the user and the doctor. Furthermore, saidcoefficient are not an average or a statistical but can easily bemeasured precisely for each diabetic user.

In one embodiment, the first term is a function of the specific insulinsensitivity factor and/or the second term is a function of both thespecific insulin response time and the specific insulin sensitivityfactor.

In one embodiment, the first term is a function of:

$\frac{{x_{1}\left( {nT}_{s} \right)} - x_{1_{ref}}}{\theta_{2}}$

In one embodiment, the second term is a function of:

θ₃×(x₂(nT_(s))+x₃(nT_(s)));

wherein:

-   -   x₂(nT_(s)) is the plasma insulin rate; and    -   x₃(nT_(s)) is the subcutaneous insulin rate.

In one embodiment, x₂(nT_(s)) is computed asx₂(nT_(s))=U_(Bas)−1/θ₂×{dot over (x)}₁(nT_(s)) or x₂(nT_(s)) is afunction of U_(Bas)−1/θ₂×{dot over (x)}₁(nT_(s)).

In one embodiment, x₃(nT_(s)) is computed asx₃(nT_(s))=x₂(nT_(s))−θ₃×({umlaut over (x)}₁(nT_(s)))/θ₂ or x₃(nT_(s))is a function of x₂(nT_(s))−θ₃×({umlaut over (x)}₁(nT_(s)))/θ₂.

{dot over (x)}₁(nT_(s)) is the time derivative of the blood glucoselevel x₁(nT_(s)) and {umlaut over (x)}₁(nT_(s)) is the second timederivative of the blood glucose level x₁(nT_(s)). U_(Bas) is an user'sspecific basal insulin injection rate.

In one embodiment, the insulin dose to be injected comprises at least aproportional component to the blood glucose level, a derivativecomponent to the blood glucose level and a second derivative componentto the blood glucose level.

In one embodiment, the insulin dose to be injected does not comprise aterm which is function of an integral of the blood glucose level.

In one embodiment, the insulin dose to be injected comprises a fourthterm being a function of a second insulin dose corresponding to the doseof insulin to be injected compensating a predefined ingested quantity ofglucose by the user. The advantage is to take into account an amount ofglucose ingested by the user during the day or during the method.

In one embodiment, the corrective factor is positive or strictlypositive and strictly inferior to 1.

In one embodiment, the step of computing an insulin dose is executed bya calculator.

In one embodiment, the method is implemented by a computer.

In one embodiment, the method further comprises the step of transmittingthe computed insulin dose to be injected to the insulin injectiondevice.

According to a eleventh aspect, the invention further relates to asystem for delivering insulin, said system comprising:

-   -   a processor comprising instructions to operate the method        according to the tenth aspect of the present invention;    -   an insulin injection device; and    -   a sensor for measuring the blood glucose level of an user.

In one embodiment, the system further comprises a transmitter totransmit data from the sensor to the processor and to transmit data fromthe processor to the insulin injection device.

In one embodiment, the system further comprises an interface configuredto define the at least one following parameter: a specific insulinresponse time; and/or a specific insulin sensitivity factor; and/or aspecific user basal insulin injection rate.

DETAILED DESCRIPTION

This invention proposes a method, a computer program and a systemimplementing a state feedback control law, derived from functionalinsulin therapy, in order to compensate high glycemia levels during afasting period or in a hybrid closed-loop. This state feedback controllaw computes basal-boluses injections, provides predictions on glucosedynamics using a long-term model, guarantees positivity of the control,and allows avoiding hypoglycemic episodes.

The system of the invention also offers the advantage that it is easy toset-up.

According to the invention, the tuning of the control law isindividualized simply using patient's own parameters such as for examplethe correction factor and the duration of insulin action. Thanks to theuse of the patient's own parameters, the tuning is readilyunderstandable to physicians, pump manufacturers, and patientsthemselves.

Insulin on Board

In this section, the model of the glucose-insulin dynamics used in thepresent invention is presented. Then it is established that the Insulinon Board can be computed as a combination of the states.

A long-term model of the glucose-insulin dynamics for type-1 diabetes isused in the present invention. Considering a fasting scenario, x₁ is theBG, x₂ and x₃ are the plasma and subcutaneous compartment insulin rate[U/min], respectively. The input u_(i) is the insulin injection rate[U/min]. θ₁ is the net balance between the endogenous glucose productionand the insulin independent consumption, θ₂ is the ISF and θ₃ is thetime constant of the insulin subsystem related to the DIA. The model is:

$\begin{matrix}{{{\overset{.}{x}}_{1} = {\theta_{1} - {\theta_{2}x_{2}}}},} & (5) \\{{{\overset{.}{x}}_{2} = {{{- \frac{1}{\theta_{3}}}x_{2}} + {\frac{1}{\theta_{3}}x_{3}}}},} & (6) \\{{\overset{.}{x}}_{3} = {{{- \frac{1}{\theta_{3}}}x_{3}} + {\frac{1}{\theta_{3}}{u_{i}.}}}} & (7)\end{matrix}$

Notice that all the states x_(i) and the control variable u_(i)represent physiological entities, therefore all are positive variables.

The insulin injection rate u_(i) is mostly the sum of a basal rate andboluses: u_(i)=U_(Bas)+u_(bol). Thus the states x₂ and x₃ can also bewritten as sums:

x ₃ =x ₃ _(Bas) +x ₃ _(bol) =x ₃ _(Bas) +{tilde over (x)} ₃,   (8)

x ₂ =x ₂ _(Bas) +x ₂ _(bol) =x ₂ _(Bas) +{tilde over (x)} ₂.   (9)

In fasting period, the correct basal insulin rate is established whenglycemia is maintained constant. The equilibrium values x_(2Bas) andx_(3Bas) are:

$\begin{matrix}{{U_{Bas} = {\left. \frac{\theta_{1}}{\theta_{2}}\Rightarrow x_{3_{Bas}} \right. = {x_{2_{Bas}} = \frac{\theta_{1}}{\theta_{2}}}}},} & (10)\end{matrix}$

and by using Eqs. (5) and (10) glycemia dynamics becomes:

$\begin{matrix}{{{\overset{.}{x}}_{1} = {{\theta_{1} - {\theta_{2}x_{2}}} = {\theta_{1} - {\theta_{2}\left( {x_{2_{Bas}} + {\overset{\sim}{x}}_{2}} \right)}}}},{{\overset{.}{x}}_{1} = {\theta_{1} - {\theta_{2}\left( {\frac{\theta_{1}}{\theta_{2}} + {\overset{\sim}{x}}_{2}} \right)}}},{{\overset{.}{x}}_{1} = {{- \theta_{2}}{{\overset{\sim}{x}}_{2}.}}}} & (11)\end{matrix}$

A physiological definition of Insulin on Board is either: the units ofinsulin from previous boluses that are still active in the body, or theamount of insulin in the subcutaneous and the plasma compartments afterboluses. According to the first definition, the state representation andthe input u_(bol), the IOB can be written as:

IOB(t)=∫₀ ^(t)(u _(bol)(τ)−{tilde over (x)} ₂(τ))dτ.   (12)

Now, merging Eqs. (6) and (7):

θ₃({tilde over ({dot over (x)})} ₃ +{tilde over ({dot over (x)})} ₂)=u_(bol) −{tilde over (x)} ₂.   (13)

Considering that no bolus was made before t=0, one gets {tilde over(x)}₃(0)=0 and {tilde over (x)}₂(0)=0. Then with (12) and (13):

IOB(t)=θ₃∫₀ ^(t)({tilde over ({dot over (x)})} ₃(τ)+{tilde over ({dotover (x)})} ₂(τ))dτ,

IOB(t)=θ₃({tilde over (x)} ₃(t)+{tilde over (x)} ₂(t)).   (14)

which agrees with the second physiological definition.

Another equivalent interpretation is:

IOB(t)=∫_(t) ^(∞) {tilde over (x)} ₂(τ)dτ,   (15)

when considering only previous boluses, the assumption is made thatu_(bol)(τ)=0 for any τ≥t. Thus

∫₀ ^(∞) u _(bol)(τ)dτ=∫ ₀ ^(∞) {umlaut over (x)} ₂(τ)dτ=∫ ₀ ^(t) u_(bol)(τ)dτ

IOB(t)=∫_(t) ^(∞) {umlaut over (x)} ₂(τ)dτ=∫ ₀ ^(t)(u _(bol)(τ)−{tildeover (x)} ₂(τ))dτ.   (16)

Integrating Eq. (11) and comparing it with (15)

θ₂IOB(t)=x ₁(t)−x _(∞).   (17)

which reads as the foreseen drop of glycemia level due to on boardinsulin, in other words, IOB provides long-term prediction on glycemia.

Control Law Design

According to this invention, a new control law called ‘Dynamic BolusCalculator’ (DBC) is introduced. In one embodiment, the DBC is based onthe correction bolus formula (4) with U_(Carb)=0 (i.e. considering nomeal),

$\begin{matrix}{{{{{Eq}.\mspace{14mu} (14)}\mspace{14mu} {and}\mspace{14mu} (17)},{{{with}\mspace{14mu} {ISF}} = {{CF} = {\theta_{2}\text{:}}}}}{{{\overset{\sim}{u}(t)} = {{U_{BG}(t)} - {{IOB}(t)}}},{\left. \Rightarrow{\overset{\sim}{u}(t)} \right. = {\frac{{x_{1}(t)} - x_{1{ref}}}{\theta_{2}} - {\theta_{3}\left( {{{\overset{\sim}{x}}_{3}(t)} + {{\overset{\sim}{x}}_{2}(t)}} \right)}}}}} & \; \\{{{Defining}\mspace{14mu} {{\overset{\sim}{x}}_{1}(t)}} = {{x_{1}(t)} - x_{1{ref}}}} & \; \\{{\overset{\sim}{u}(t)} = {{{\overset{\sim}{u}}_{1}(t)} = {\frac{{\overset{\sim}{x}}_{1}(t)}{\theta_{2}} - {\theta_{3}{{\overset{\sim}{x}}_{2}(t)}} - {\theta_{3}{{{\overset{\sim}{x}}_{3}(t)}.}}}}} & (18)\end{matrix}$

In one embodiment, the invention consists to use the equation (18) incontinuous. The θ₂ and θ₃ parameters are provided to the computerprogram and are tools usually handled by the patient. In consequence, anadvantage is the method according to the present invention ispersonalized and very simple to be applied to different diabetic users.

According to one embodiment, the computed global insulin injection ratecomprises a constant insulin injection rate such as basal rate and avariable insulin injection rate.

The global injection rate u_(i)(t) will be the state feedback modulatingũ₁(t) the constant insulin injection rate U_(Bas):

u _(i)(t)=U _(Bas) +ũ ₁(t).   (19)

Thus, with (6), (7), (11) and (18), the following closed loop will bestudied:

{tilde over ({dot over (x)})}(t)=A{umlaut over (x)}(t)+Bũ _(k)(t),{umlaut over (x)}(0)={dot over (x)} ₀.   (20)

ũ _(k)(t)=F _(k) {tilde over (x)},   (21)

where ũ_(k) is defined as ũ_(k)=kũ. This feedback defines an entirefamily of DBC controllers, which are part of this invention. Thematrices A, B, and F_(k) are

${A = \begin{pmatrix}0 & {- \theta_{2}} & 0 \\0 & {- \frac{1}{\theta_{3}}} & \frac{1}{\theta_{3}} \\0 & 0 & {- \frac{1}{\theta_{3}}}\end{pmatrix}},\mspace{14mu} {B = \begin{pmatrix}0 \\0 \\\frac{1}{\theta_{3}}\end{pmatrix}},\mspace{14mu} {F_{k} = {{k\begin{pmatrix}\frac{1}{\theta_{2}} & {- \theta_{3}} & {- \theta_{3}}\end{pmatrix}}.}}$

An interesting property of this family of controllers is that the totalquantity of injected insulin does not depend on k:

∫₀ ^(∞) ũ _(k)(τ)dτ=∫ ₀ ^(∞) ũ ₁(τ)dτ=ũ ₁(0).   (22)

This allows to stretch the input trajectory and to keep constant thetotal quantity of insulin injected, and k just tunes the velocity ofdecrease of glycemia without falling into hypoglycemia levels, as shownhereafter.

Input/State Positivity

In this section, important properties as stability and positivity of theclosed-loop trajectories are addressed. It is proven that this feedbackgenerates a positive control, which ensures the positivity of, {tildeover (x)}₁ that is x₁(t)≥x_(1ref).

In medical terms, this property is a guaranty of no hypoglycemicepisodes.

According to Eqs. (20-21), the closed-loop system reads as:

$\begin{matrix}{{{\overset{.}{\overset{\sim}{x}}(t)} = {{\left( {A + {BF}_{k}} \right){\overset{\sim}{x}(t)}} = {\overset{\sim}{A}{\overset{\sim}{x}(t)}}}},{{\overset{.}{\overset{\sim}{x}}(t)} = {\begin{pmatrix}0 & {- \theta_{2}} & 0 \\0 & {- \frac{1}{\theta_{3}}} & \frac{1}{\theta_{3}} \\\frac{k}{\theta_{2}\theta_{3}} & {- k} & {{- k} - \frac{1}{\theta_{3}}}\end{pmatrix}{\overset{\sim}{x}(t)}}},} & (23)\end{matrix}$

which is a stable system with eigenvalues

$\lambda_{1} = {\lambda_{2} = {- \frac{1}{\theta_{3}}}}$

and λ₃=−k, for all k>0.

The positivity of the input/state trajectories, i.e. ũ_(k)(t)≥0 and{tilde over (x)}(t)≥0 ∀t≥0, is discussed through the notion ofpositively invariant sets.

-   -   Definition 1: Given a dynamical system {dot over (x)}−Dx. x∈        ^(n), and a trajectory x(t, x₀), where x₀ is a initial        condition, a nonempty set M ⊆        ^(n) is a positive invariant set if ∀x₀ ∈ M ⇒ x(t, x₀)∈M ∀t≥0.    -   Definition 2: If G∈        ^(r×n) then M(G) denotes the polyhedron

M(G)={x∈

^(n) |Gx≥0}.   (24)

Proposition 1: The polyhedral set M(G) is a positively invariant set forthe system of Definition 1 if and only if there exists a Metzler matrix

H∈

^(r×r), i.e. H_(ij)≥0 for i≠j, such that:

GD−HG=0.   (25)

Note that the polyhedron defined by the positive orthant in R³ is not apositively invariant set for the system (23) since in this case G=I andD=Ã, and H must be Ã; however the latter matrix is not Metzler.

In order to find the maximal invariant set, the system (23) istransformed to its Jordan form, i.e.

ż(t)=P ⁻¹ ÃPz(t)=Jz(t).   (26)

where {tilde over (x)}=Pz, and

$\begin{matrix}{J = \begin{pmatrix}{- \frac{1}{\theta_{3}}} & 1 & 0 \\0 & {- \frac{1}{\theta_{3}}} & 0 \\0 & 0 & {- k}\end{pmatrix}} & (27)\end{matrix}$

through the matrix of change of coordinates

$\begin{matrix}{P = {\begin{pmatrix}{\theta_{2}\theta_{3}} & {\theta_{2}\theta_{3}^{2}} & 1 \\1 & 0 & \frac{k}{\theta_{2}} \\0 & \theta_{3} & \frac{k\left( {1 - {k\; \theta_{3}}} \right)}{\theta_{2}}\end{pmatrix}.}} & (28)\end{matrix}$

Notice that the positive orthant in the new coordinates is a positivelyinvariant set because the matrix J is Metzler.

The state trajectories {tilde over (x)}(t) in z-coordinates are

$\begin{matrix}{{{\overset{\sim}{x}}_{1}(t)} = {{\theta_{2}\theta_{3}{e^{\frac{- t}{\theta_{3}}}\left( {{z_{1}(0)} + {\left( {\theta_{3} + t} \right){z_{2}(0)}}} \right)}} + {e^{- {kt}}{z_{3}(0)}}}} & (29) \\{{{{\overset{\sim}{x}}_{2}(t)} = {{e^{\frac{- t}{\theta_{3}}}\left( {{z_{1}(0)} + {{tz}_{2}(0)}} \right)} + {e^{- {kt}}\frac{k}{\theta_{2}}{z_{3}(0)}}}},} & (30) \\{{{{\overset{\sim}{x}}_{3}(t)} = {{e^{\frac{- t}{\theta_{3}}}\theta_{3}{z_{2}(0)}} - {e^{- {kt}}\frac{k\left( {{k\; \theta_{2}} - 1} \right)}{\theta_{2}}{z_{3}(0)}}}},} & (31)\end{matrix}$

and the control trajectory becomes

$\begin{matrix}{{{{\overset{\sim}{u}}_{k}(t)} = {e^{- {kt}}\frac{{k\left( {{k\; \theta_{3}} - 1} \right)}^{2}}{\theta_{2}}{z_{3}(0)}}},} & (32)\end{matrix}$

It verifies that

$\begin{matrix}{{\frac{{k\left( {{k\; \theta_{3}} - 1} \right)}^{2}}{\theta_{2}}z_{3}} = {{\overset{\sim}{u}}_{k} = {{k\left( {{\frac{1}{\theta_{2}}{\overset{\sim}{x}}_{1}} - {\theta_{3}\left( {{\overset{\sim}{x}}_{2} + {\overset{\sim}{x}}_{3}} \right)}} \right)}.}}} & (33)\end{matrix}$

Now, the property (22) can be proved. Consider the integral of Eq. (32)

${\int_{0}^{\infty}{{{\overset{\sim}{u}}_{k}(\tau)}d\; \tau}} = {{{{- \frac{\left( {{k\; \theta_{3}} - 1} \right)^{2}}{\theta_{2}}}e^{- {kt}}{z_{3}(0)}}|_{0}^{\infty}} = {\frac{\left( {{k\; \theta_{3}} - 1} \right)^{2}{z_{3}(0)}}{\theta_{2}}.}}$

Replacing z₃ from Eq. (33) in the latter equation

${\int_{0}^{\infty}{{{\overset{\sim}{u}}_{k}(\tau)}d\; \tau}} = {{{\frac{1}{\theta_{2}}{{\overset{\sim}{x}}_{1}(0)}} - {\theta_{3}\left( {{{\overset{\sim}{x}}_{2}(0)} + {{\overset{\sim}{x}}_{3}(0)}} \right)}} = {{{\overset{\sim}{u}}_{1}(0)}.}}$

As the control trajectory Eq. (32) is an exponential function dependingon k, that allows us to stretch the trajectory ensuring that the samequantity of insulin is administered for all k>0.

The following theorem restates the positivity of the first orthant inz-space, but in the x-coordinates.

Theorem 1: Consider the sets M₁=

{tilde over (x)}∈

³|{tilde over (x)}≥0}, and

$M_{2} = {\left\{ {\overset{\sim}{x} \in {\mathbb{R}}^{3}} \middle| {{k\left( {{\frac{1}{\theta_{2}}{\overset{\sim}{x}}_{1}} - {\theta_{3}\left( {{\overset{\sim}{x}}_{2} + {\overset{\sim}{x}}_{3}} \right)}} \right)}^{\prime} \geq 0} \right\}.}$

The maximal positively invariant polyhedron of system (23) is

M=M₁ ∩ M₂.   (34)

Proof: It is clear that the condition x(0)≥0 is necessary but it doesnot ensure that the {tilde over (x)}-trajectories remain positive fort>0, because the matrix Ã is not Metzler. However, the set {{tilde over(x)}∈

³|{tilde over (x)}≥0} contains any positively invariant set, i.e.

M ⊂ M₁.

In z-coordinates, the positively invariant set is the first orthant in

³, and by Eq. (32), z₃ is proportional to ũ, therefore ũ≥0 is anecessary condition but not sufficient. Then

M ⊂ M₂.

Now, using Definition 2, the polyhedron {Gx≥0} where

$\begin{matrix}{{G = \begin{pmatrix}1 & 0 & 0 \\0 & 1 & 0 \\0 & 0 & 1 \\\frac{k}{\theta_{2}} & {{- k}\; \theta_{3}} & {{- k}\; \theta_{3}}\end{pmatrix}},} & (35)\end{matrix}$

characterizes the positive invariant set for the system (23), andtranslates the positive orthant in z to x-coordinates. That is verifiedusing Proposition 1 with the following H-matrix

$\begin{matrix}{{H = \begin{pmatrix}{- \frac{1}{\theta_{3}}} & 0 & \theta_{2} & \frac{\theta_{2}}{k\; \theta_{3}} \\0 & {- \frac{1}{\theta_{3}}} & \frac{1}{\theta_{3}} & 0 \\0 & 0 & {- \frac{1}{\theta_{3}}} & \frac{1}{\theta_{3}} \\0 & 0 & 0 & {- k}\end{pmatrix}},} & (36)\end{matrix}$

And D=Ã. As this Proposition is a necessary and sufficient condition,then the set M=M₁ ∩ M₂ is the maximal positively invariant polyhedronfor system (23).

According to Definition 1, the nonempty set M ⊆

³ is the positively invariant polyhedron of the system (22) controlledby Eq. (21), that is, if the system starts inside M, it will remainsthere for any t>0. As the insulin subsystem is indeed positive, thecondition to ensure the positivity can be summarized as ũ≥0.

From a medical point of view, the positivity of the input ensures that

≥0, i.e. guaranties the exclusion of hypoglycemia episodes:

(x ₁ ≥x ₁ _(ref) , ∀t≥0):

Moreover, positivity of the control stands in agreement with themanagement of insulin injection.

Nonetheless, the eigenvalues of the insulinemia subsystem

$\left( {- \frac{1}{\theta_{3}}} \right)$

are not modified by the control law. Consequently the performance of theclosed-loop depends on the patient's θ₃ parameter.

Robustness

Robustness is a decisive issue as it ensures that the controller willwork safely on the non-nominal diabetic patient.

According to one embodiment, the processor for computing further definesa reference level of Blood Glucose; and wherein at the final endpoint ofeach time interval, the global insulin injection rate is corrected,taking into account the gap between the measured level of Blood Glucoseand a reference level of Blood Glucose.

Delay Margin

Delays are well known to destabilize closed-loops. Here, it is studiedthe robustness with respect to delays that were not taken into accountin the model. These delays appear naturally in the closed-loop of theartificial pancreas. Analyzing the state feedback (23), the target looptransfer is given by:

$\begin{matrix}{{L_{Target} = {F_{k} \cdot \left( {{sI} - A} \right)^{- 1} \cdot B}}{L_{Target} = \frac{k}{s}}} & (37)\end{matrix}$

The phase margin of L_(Target) is

$\frac{\pi}{2}$

at pulsation ω_(k)=k. This leads to a delay margin

${M_{r} = \frac{\pi}{2\omega_{k}}},{M_{r} = \frac{\pi}{2k}},$

which is the maximum added delay that does not destabilize the loop. Forinstance, setting Mr=25 min resolves

k=k _(r)≅10⁻³ rad/s.

Performance by Regulation

The closed-loop shows that it is possible to reject disturbance actingas an output step, but not as a ramp. In the case of a ramp disturbance,the speed error will be:

$\epsilon = \frac{slope}{k}$

Parameters Uncertainties

Because the behavior of the organism of the user could be different ofthe nominal behavior, the calculated bolus is not delivered in one dosebut is spread in time. In this case, the steep fall in Blood Glucoserate is limited.

With the factor k_(r) which ensures stability of the closed-loop fordelays lower than 25 minutes, robustness with respect to parametersuncertainties is studied. The state feedback is:

$\begin{matrix}{{{\hat{F}}_{k_{r}} = {k_{r}\left\lbrack {\frac{1}{{\hat{\theta}}_{2}} - {\hat{\theta}}_{3} - {\hat{\theta}}_{3}} \right\rbrack}},} & (38)\end{matrix}$

Where {circumflex over (θ)}_(l) are estimates of the model parametersθ_(i). Considering the state feedback A+B{circumflex over (F)}_(k) _(r), the target loop transfer is given by:

$\begin{matrix}{{\hat{L} = {{{\hat{F}}_{k_{r}}\left( {{sI} - A} \right)}^{- 1} \cdot B}}{\hat{L} = {k_{r}\frac{\frac{1}{\theta_{3}^{3}}\left( {\frac{\theta_{2}}{\theta_{2}} + {s\; {\hat{\theta}}_{3}} + {\theta_{3}{\hat{\theta}}_{3}{s\left( {s + \frac{1}{\theta_{3}}} \right)}}} \right)}{{s\left( {s + \frac{1}{\theta_{3}}} \right)}^{2}}}}} & (39)\end{matrix}$

With the Nyquist criterion, FIG. 1 shows that stability is ensured evenwith great parameters uncertainties. Moreover, the delay margin is stillgood as it is equal to 12 min at worst for {circumflex over(θ)}_(2,3)=50%×θ_(2,3).

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 represents L_(Target) and {circumflex over (L)} with50%×θ_(2,3)≤{circumflex over (θ)}_(2,3)≤200%×θ_(2,3).

FIG. 2 represents state feedback F_(kr) with delay T_(r) not taken intoaccount and with well-known model parameter ({circumflex over (θ)}=θ).

FIG. 3 represents Dawn Phenomenon: open and closed-loop with statefeedback F_(kr).

FIG. 4 represents closed-loop with state feedback F_(kr) and CFunderestimated.

FIG. 5 represents Dynamic Bolus Calculator control law with Adult01 fromUVA/Padova T1DMS wearing Generic pump and Dexcom 70 CGM device.Estimated parameters: {circumflex over (θ)}₁=1.6 mg/dl/min, {circumflexover (θ)}₂=75 mg/dl/U, {circumflex over (θ)}₃=38 min.

FIG. 6 represents the glycemia level of a diabetic patient (above) andthe amount of insulin injection (below) during time. Two models arerepresented: the model with one bolus of insulin injection (named“bolus” on the graph) and the model according to one embodiment of thepresent invention wherein the insulin injection is continuouslycalculated and injected according to the present invention and whereinthe k coefficient is strictly inferior to 1 (“spread bolus”) and whereinthe loop is closed at t=30 min.

FIG. 7 represents an enlargement of the graph on FIG. 6.

RESULTS

The following simulations are conducted under meal-free scenarios. Thereference is set to 100 mg/dl and loop is closed at t=60 minutes.

The patient's parameters, derived from a real patient, are θ₁=0.85mg/dl/min, θ₂=70 mg/dl/U and θ₃=62 min.

Closed-Loop with Delay

All CGM devices introduce some delay due to the physio-logical time lagbetween blood glucose and interstitial glucose concentration.

FIG. 2 illustrates the closed-loop (23) with a delay T_(r) added to thestate {tilde over (x)}₁. The state feedback F_(kr) uses the delayedoutput {tilde over (x)}₁(t−T_(r)) and the current states {tilde over(x)}₂(t) and {tilde over (x)}₃(t).

At the beginning (t=0) BG=300 mg/dl, IOB=0 U.

As argued in section “Input/State Positivity”, in the nominal case(Tr=0):

-   -   the control remains positive ũ_(kr)≥0 (as u_(i)≥U_(Bas));    -   the state remains positive {tilde over (x)}₁≥0 (as x₁≥x_(ref)),        there is no hypoglycemic episode;    -   the speed of return in euglycemia depends on θ₃ as explained in        section “Input/State Positivity”.

When T_(r)=15 min

-   -   loop remains stable as expected (see part “Delay Margin”);    -   there is no hypoglycemic episode;    -   although it was no demonstrated in case of delay ∫ũ_(k) _(r)        |_(τ) _(r) ₌₀=∫ü_(k) _(r) |_(τ) _(r) ₌₁₅=2.86 U which        corresponds to the injection advised by a bolus wizard at t=1 h        with respect to Eq. (4).

The Dawn Phenomenon

The dawn phenomenon is an increase in blood glucose in the night due tosurge in hormone secretion. FIG. 3 illustrates the closed-loop (23)using the state feedback F_(kr), with the dawn phenomenon modeled by aramp disturbance on the output of 25 mg/dl/h between 2 and 6 a.m.

In open loop, under constant basal rate, the dawn phenomenon bringsglycemia to 200 mg/dl at 8 a.m.

In closed-loop:

-   -   the state feedback control law produces a temporary basal        delivery rate;    -   as expected (see part “Performance by regulation”) the effect of        the ramp disturbance;    -   and BG recovers euglycemia at 8 a.m.

Parameter Uncertainties

In this section uncertainty on CF will be addressed. Assuming that thepatient has CF of 70 mg/dl/U, the worst case is considered: CF isunderestimated ({circumflex over (θ)}₂=70%×θ₂=50 mg/dl/U).

The initial glycemia is 300 mg/dl and the target is 100 mg/dl. In openloop, this would involve dramatic consequences as the computed bolus((300−100)/50=4 U) should lower the glycemia by CF×4 U=280 mg/dl andlead the patient to severe hypoglycemia (BG above 20 mg/dl).

FIG. 4 shows high safety of the closed-loop as despite anunderestimation of the CF, the glycemia reaches target with nohypoglycemia (the minimum BG is 96 mg/dl).

UVA/Padova T1DM Simulator

The distributed version of the UVA/Padova has been approved by the Foodand Drug Administration as a preclinical testing platform for controlalgorithm. With several virtual patients, it also includes models ofpumps and CGM devices. This simulator is used to demonstrate the safetyand efficiency of the DBC control algorithm. The parameters ({circumflexover (θ)}₁, {circumflex over (θ)}₂ and {circumflex over (θ)}₃) of thevirtual patient are identified from a previous scenario. The initial BGis 300 mg/dl. At t=0 the loop is closed. The virtual patient uses aGeneric pump (increase step is 0.05 U) and a CGM device (Dexcom 70)which introduces delay and noise and has a sampling time of 5 minutes.The reference has been set to 110 mg/dl.

FIG. 5 shows a good performance of the closed-loop as:

-   -   BG is into euglycemia at 2 h 3;    -   there is no hypoglycemia (minimum BG is 82 mg/dl).

Robustness of the closed-loop is also demonstrated as:

-   -   the model of the virtual patient is non-nominal;    -   the CGM devices introduces delays, noise and a sampling time of        5 min;    -   the pump has a minimum delivery rate step of 0.05 U;    -   there was no saturation of the controller.

CONCLUSION

Individualization of the controller and accurate prediction for MPCalgorithm are still an open problem in the artificial pancreas project.Here, a novel closed-loop is developed. With a structure derived frombolus advisors, this controller is simply tuned with individualizedcharacteristics of the patient (CF and DIA). Thus it is immediatelycomprehensive to physicians, and patients.

The main feature of this control law is to ensure the positivity oftrajectories. This guarantees that the glycemia remains greater than itsreference, at least in the nominal case and allows the controller tocope with the positivity constraint of the insulin injection.

As in practice there are some delays, and parameter uncertainties, arobustness analysis is added.

Finally, through simulations the performance of the loop is assessed,for the nominal case, and for a more realistic scenario, the UVA/Padovasimulator was used to implement this.

As the dynamics of the insulinemia subsystem are not modified, a controllaw is envisioned to accelerate the response. Also, it is necessary totackle the problem with meals and a long term scenario. The goodperformance obtained with the simulator encourages to propose clinicaltrials in this subject.

EXAMPLE 1 Comparison Between One Single Bolus Dose and a Portioned BolusDose

One degree of freedom in the tuning of the controller is the time inwhich a bolus is delivered. One can choose to deliver a bolus in onesingle dose or to partition said single dose. The last point is a pledgeof security. FIG. 6 shows a simulation wherein patient parameters areknown. The loop is closed after a duration of 30 minutes. One willnotice that a “Bolus” instruction generates one single bolus dose thendelivers the basal bolus dose (basal rate). The “Spread Bolus” workswith a control period of 15 minutes and delivers 77% of the bolus in 1hour.

FIG. 6 shows a simulation when the patient compensatory value is notwell entered. The patient, here, has a compensatory value equal to 70mg/dl/U and the controller works with a wrong value of 50 mg/dl/U. Thus,the generated bolus dose is of 4 U, which leads the glycemia to a finalvalue of 300−4*70 that corresponds to a calculated value of 20 mg/dl.

However, the controller observes the deviation between the measuredglycemia value and the target glycemia value, given by the IOB andcorrects at the next instruction by retracting a part of the basal dose(1.1 U for 4 hours). The global value remains strictly positive. Asmentioned, the “spread bolus” instruction provides a pledge of security.Indeed, 3.5 U are injected in 1 h 45 then 0.6 U are retracted of thebasal dose.

One can observe on the enlargement on FIG. 7, with the “bolus”instruction that the corrector response to the overdose is performed atthe second clock tick in order to reach its maximal effect (75% of theretracted basal dose) one hour after the bolus release. Therefore, the“spread bolus”, expressed by a coefficient k_(d) strictly inferior to 1,allows keeping better margins (45% of the basal value retracted in theworst case).

1. A method for controlling an insulin injection device of a user,comprising iteratively the steps of: determining a time interval T_(s);receiving the blood glucose level x₁(nT_(s)); computing an insulin doseto be injected u(nT_(s)) in the next time interval; wherein the insulindose to be injected u(nT_(s)) comprises at least: a first term being afunction of the comparison between the received blood glucose levelx₁(nT_(s)) and a predefined blood glucose level target x₁ _(ref) in apreliminary step of the method; a second term, the second term being anestimated value of the insulin dose still active in the body IOB(nT_(s))of the user; the second term being superior or equal to the first termat each iteration of the method.
 2. The method according to claim 1,wherein the first and the second terms of the insulin dose to beinjected u(nT_(s)) is a function of a corrective factor (k_(d)) inferioror equal to 1, the corrective factor (k_(d)) being configured to adaptthe duration of the injection to a predefined duration reference.
 3. Themethod according to claim 2, wherein the first and the second terms ofthe insulin dose to be injected u(nT_(s)) is a linear function of thecorrective factor (k_(d)).
 4. The method according to claim 1, whereinthe insulin dose to be injected u(nT_(s)) comprises: a third termcalculated on at least one specific injection rate of a predefined userprofile.
 5. The method according to claim 1, wherein the insulin dose tobe injected u(nT_(s)) is a function of at least one of the followingpredefined user profile parameters: a specific insulin response time θ₃;and/or a specific insulin sensitivity factor θ₂.
 6. The method accordingto claim 1, wherein the first term is a function of:$\frac{{x_{1}({nTs})} - x_{1_{ref}}}{\theta_{2}}$ wherein θ₂ is aspecific insulin sensitivity factor.
 7. The method according to claim 1,wherein the second term is a function of:θ₃×(x₂(nT_(s))+x₃(nT_(s))) wherein θ₃ is a specific insulin responsetime; x₂(nT_(s)) is the plasma insulin rate; and x₃(nT_(s)) is thesubcutaneous insulin rate.
 8. The method according to claim 7, whereinx₂(nT_(s)) is computed as x₂(nT_(s))=U_(Bas)−1/θ₂×{dot over(x)}₁(nT_(s)); and x₃(nT_(s)) is computed asx₃(nT_(s))=x₂(nT_(s))−θ₃×({umlaut over (x)}₁(nT_(s)))/θ₂; wherein {dotover (x)}₁(nT_(s)) is the time derivative of the blood glucose levelx₁(nT_(s)) and {umlaut over (x)}₁(nT_(s)) is the second time derivativeof the blood glucose level x₁(nT_(s)), and further wherein U_(Bas) is auser's specific basal insulin injection rate and θ₂ is a specificinsulin sensitivity factor.
 9. The method according to claim 1, whereinthe insulin dose to be injected u(nT_(s)) comprises at least aproportional component to the blood glucose level x₁(nT_(s)), aderivative component to the blood glucose level x₁(nT_(s)) and a secondderivative component to the blood glucose level x₁(nT_(s)).
 10. Themethod according to claim 1, wherein the insulin dose to be injectedu(nT_(s)) does not comprise a term which is a function of an integral ofthe blood glucose level x₁(nT_(s)).
 11. The method according to claim 1,wherein the insulin dose to be injected u(nT_(s)) comprises: a fourthterm being a function of a second insulin dose u_(Carb) corresponding tothe dose of insulin to be injected compensating a predefined ingestedquantity of glucose by the user.
 12. The method according to claim 2,wherein the corrective factor (k_(d)) is positive and strictly inferiorto
 1. 13-14. (canceled)
 15. The method according to claim 1, said methodfurther comprising the step of: transmitting the computed insulin doseto be injected u(nT_(s)) to the insulin injection device. 16-17.(canceled)
 18. A computer-implemented method for controlling an insulininjection device of a user, comprising iteratively the steps of:determining a time interval T_(s); receiving the blood glucose levelx₁(nT_(s)); computing an insulin dose to be injected u(nT_(s)) in thenext time interval; wherein u(nT_(s))=k_(d)×ũ(nT_(s))+T_(s)×U_(Bas)wherein: k_(d) is a tuning parameter strictly positive and inferior orequal to 1; ũ(nT_(s)) is a correction insulin dose computed asũ(nT_(s))=u_(BG)(nT_(s))−IOB(nT_(s)); wherein u_(BG)(nT_(s)) is theinsulin dose needed to reach the blood glucose level x₁(nT_(s)) to ablood glucose level target x₁ _(ref) without considering previousinsulin injections; and IOB(nT_(s)) is the insulin dose still active inthe body of the user; and U_(Bas) is the user's specific basal insulininjection rate. transmitting the computed insulin dose to be injectedu(nT_(s)) to the insulin injection device. 19-38. (canceled)
 39. Acomputer program comprising instructions for the implementation of stepsof a method according to claim 1 when the program is executed by acomputer.
 40. A system for delivering insulin, the system comprising: aprocessor comprising instructions to operate the method according toclaim 1; an insulin injection device; and a sensor for measuring theblood glucose level of a user.
 41. The system according to claim 40further comprising an interface configured to define the at least onefollowing parameter: a specific insulin response time θ₃; and/or aspecific insulin sensitivity factor θ₂; and/or a specific user basalinsulin injection rate U_(Bas).
 42. A system for delivering insulin, thesystem comprising: an insulin injection device; a sensor for measuringthe blood glucose level of a user; a processor comprising instructionsconfigured to iteratively operate steps of: determining a time intervalT_(s); receiving the blood glucose level x₁(nT_(s)); computing aninsulin dose to be injected u(nT_(s)) in the next time interval; whereinthe insulin dose to be injected u(nT_(s)) comprises at least: a firstterm being a function of the comparison between the received bloodglucose level x₁(nT_(s)) and a predefined blood glucose level target x₁_(ref) in a preliminary step of the method; a second term, the secondterm being an estimated value of the insulin dose still active in thebody IOB(nT_(s)) of the user; the second term being superior or equal tothe first term at each iteration.
 43. The system according to claim 42,wherein the first and the second terms of the insulin dose to beinjected u(nT_(s)) is a function of a corrective factor (k_(d)) inferioror equal to 1, the corrective factor (k_(d)) being configured to adaptthe duration of the injection to a predefined duration reference. 44.The system according to claim 42, wherein the insulin dose to beinjected u(nT_(s)) is a function of at least one of the followingpredefined user profile parameters: a specific insulin response time θ₃;and/or a specific insulin sensitivity factor θ₂.